rm(list = ls(all = TRUE))
gc()
##           used (Mb) gc trigger (Mb) max used (Mb)
## Ncells  617808 33.0    1414808 75.6   686457 36.7
## Vcells 1163527  8.9    8388608 64.0  1877160 14.4

1 packages

# %>% 
library(magrittr)
library(ggplot2)


`%!in%` = Negate(`%in%`)

2 palette

n = 4
col = RColorBrewer::brewer.pal(n, 'Paired')

pal = col[seq(1, n, 2)]

3 f-ctions

3.1 summary stat

summary_stats <- function(df, measurevar, groupvars, conf.level = 0.95) {
  df %>%
    dplyr::group_by(across(all_of(groupvars))) %>%
    dplyr::summarise(
      N = sum(!is.na(.data[[measurevar]])),
      mean = mean(.data[[measurevar]], na.rm = TRUE),
      sd = sd(.data[[measurevar]], na.rm = TRUE),
      .groups = "drop"
    ) %>%
    dplyr::mutate(
      se = sd / sqrt(N),
      ci = se * qt(conf.level/2 + 0.5, N - 1)
    )
}

3.2 scale to the average of the reference group

process_data <- function(df, groupvars, measurevar, scale_treatment = "noninoculated") {

  df.long = df %>%
    tidyr::pivot_longer(cols = 3:ncol(df), names_to = "transcript",
                        values_to = measurevar, values_drop_na = TRUE)
  data.SE = summary_stats(df.long, measurevar, groupvars)
  scale_reference = data.SE %>%
    dplyr::filter(Treatment == scale_treatment) %>%
    dplyr::select(Tissue, transcript, mean) %>%
    dplyr::rename(scale_mean = mean)
  df.scaled <- dplyr::left_join(df.long, scale_reference, by = c("Tissue", "transcript")) %>%
    dplyr::mutate(scaled = .data[[measurevar]] / scale_mean) %>%
    dplyr::select(-scale_mean)
  df.scaled = df.scaled %>% 
    dplyr::arrange(dplyr::desc(transcript), dplyr::desc(Treatment))
  shoots = df.scaled %>% dplyr::filter(Tissue == "shoots")
  roots = df.scaled %>% dplyr::filter(Tissue == "roots")
  
  return(list(
    shoots = shoots,
    roots = roots
  ))
}

3.3 permutation t-test

  • for non-i.i.d, hypotheses
  • a non-parametric robust alternative to traditional t-tests
  • small sample sizes
  • robust to unequal variances
perm_test_by_transcript <- function(mydata.long, measurevar = "measurement", groupvar = "Treatment") {
  
  set.seed(123456)
  
  temp = data.frame()
  transcript_levels = levels(mydata.long$transcript)
  treatment_levels = levels(mydata.long[[groupvar]])
  treatment_pairs = combn(treatment_levels, 2, simplify = FALSE)
  
  for(i in transcript_levels) {
    data = mydata.long[mydata.long$transcript == i, ]
    k = 12 # nrow(data)
    pvalues = purrr::map_dbl(treatment_pairs, function(x) {
      subset_data = base::subset(data, data[[groupvar]] %in% x)
      res = MKinfer::perm.t.test(
        formula = stats::as.formula(base::paste(measurevar, "~", groupvar)),
        data = subset_data,
        alternative = "two.sided",
        mu = 0,
        paired = FALSE,
        var.equal = FALSE,
        conf.level = 0.95,
        perm.conf.int = 0.95,
        symmetric = TRUE,
        p.adjust.method = "BH",
        detailed = TRUE,
        # split k observations into two nontrivial groups
        R = sum(choose(k, 1:(k-1))), # alternative (k, k/2)
        set.seed = 123456
      )
      res$perm.p.value
    })

    
    tmp = base::as.data.frame(base::t(pvalues))
    colnames(tmp) = purrr::map_chr(treatment_pairs, ~base::paste(.x, collapse = ' vs '))
    rownames(tmp) = i
    
    temp <- base::rbind(temp, tmp)
  }
  
  
    return(temp)
}

3.4 plot

plot_gene_expression <- function(data.SE
                                 , data.long
                                 , stat.test.sig
                                 , transcripts_excl = c('13-LOX', 'PTI5')
                                 , color_values
                                 , facet_cols = 6
                                 , y_scales = NULL
                                 , dodge_width = 0.8,
                                 plot_title = ""
                                 ) {
  
  dodge = position_dodge(width = dodge_width)
  
  data.SE.filtered = dplyr::filter(data.SE, !transcript %in% transcripts_excl)
  data.long.filtered = dplyr::filter(data.long, !transcript %in% transcripts_excl)
  
  stat.test.filtered = NULL
  if (nrow(stat.test.sig) > 0) stat.test.filtered = dplyr::filter(stat.test.sig, !transcript %in% transcripts_excl)
  
  p = ggplot(data.SE.filtered, aes(x = Treatment, y = mean)) +
    geom_point(size=3.5, shape = 22, position = dodge, aes(fill = Treatment), colour = "black") +
    geom_point(data = data.long.filtered, aes(x = Treatment, y = scaled, fill = Treatment), colour = "black",
               size = 2.0, shape = 21,
               position = dodge) +
    geom_errorbar(aes(ymin = mean - se, ymax = mean + se), width = 0.3, lwd = 0.5,
                  position = dodge) +
    facet_wrap(~ transcript, ncol = facet_cols, scales = "free", drop = TRUE)
  
  if (!is.null(y_scales)) {
    p = p + ggh4x::facetted_pos_scales(y = y_scales)
  }
  
  p = p +
    # geom_hline(yintercept = 2.5, alpha = 0.0) +
    geom_hline(yintercept = 0, alpha = 0.0) +
    geom_hline(yintercept = 1, alpha = 0.5, linetype = "dotted", col = "gray45") +
    ggtitle(plot_title) +
    theme_bw() +
    scale_colour_manual(name = "", values = rev(color_values)) +
    scale_fill_manual(name = "", values = rev(color_values)) +
    labs(x = "", y = "Relative gene expression (+/- SE)"
         # , subtitle = "Permutation t-test measurements"
         ) +
    theme(plot.subtitle = element_text(size = 10),
          axis.text = element_text(size = 12.5),
          axis.text.x = element_text(size = 12.5, angle = 90),
          axis.title = element_text(size = 12.5, face = "bold"),
          strip.text = element_text(size = 12.5),
          title = element_text(size = 15, face = "bold"),
          # axis.ticks.x = element_blank(),
          legend.key.height = unit(1.5, "cm"),
          legend.key.width = unit(1.75, "cm"),
          legend.text = element_text(size = 12.5),
          legend.title = element_text(size = 12.5),
          legend.background = element_rect(fill = "transparent", size = 0.5, linetype = "dotted"),
          legend.position = "top",
          legend.justification = "right",
          plot.background = element_blank(),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          axis.title.y.right = element_blank(),
          axis.text.y.right = element_blank(),
          # axis.ticks.y = element_blank(),
          axis.text.y = element_text(size = 12.5, margin = margin(r = 0)),
          panel.spacing = unit(1, "lines"),
          strip.background = element_rect(size = 0.5, fill = "transparent", color = NA) ,
          panel.border = element_blank(),
          axis.line = element_line(color = "black")
          )
  
    p = p + theme(legend.position = "none")

  


  
  if (!is.null(stat.test.filtered) && nrow(stat.test.filtered) > 0) {
    
    p = p + ggpubr::stat_pvalue_manual(
      stat.test.filtered,
      label = "p.adj.signif",
      xmin = "xmin",
      xmax = "xmax",
      y.position = "y.position",
      hide.ns = FALSE,
      tip.length = 0.01,
      step.increase = 0.0,
      inherit.aes = FALSE 
    )
    
  }
    
  return(p)
  
}

3.5 distribution, variance and effects

  • on qPCR data
  • Denote: LOQ values skew distributions
  • Denote: distribution tests: low power for small sample size, see QQ of residuals
  • if most or all points fall inside the shaded confidence band, the sample’s distribution does not show strong evidence of departure from normality at the plotted sample size and confidence level
  • a few isolated points outside the band at the extremes are common with small samples and do not necessarily indicate a severe problem
  • Denote: heteroscedasticity: some tests are meant to be used with normally distributed data, but can tolerate relatively low deviation from normality; Levene’s test with mean, more robust test Brown-Forsythe with median, Fligner-Killeen when data are non-normally distributed or when problems related to outliers in the dataset cannot be resolved
  • Denote: effect size: Cohen’s d - parametric with assumptions, Wilcoxon Effect Size - the non-parametric alternative; different sensitivity to outliers
dens_and_effect <- function(mydata.long, p_colors) {
  
  # Plot density with x="measurement"
  p1 = ggpubr::ggdensity(mydata.long, x = "measurement",
              add = "mean", rug = TRUE,
              color = "Treatment", fill = "Treatment",
              facet.by = 'transcript') +
    scale_fill_manual(name = "", values = rev(p_colors)) +
    scale_color_manual(name = "", values = rev(p_colors)) +
    ggtitle("measurement") +
    facet_wrap(~transcript, ncol = 4, scales = "free")
  print(p1)
  
  # optional: Plot density with x="scaled"- sam e shape wih a shift
  # p2 = ggpubr::ggdensity(mydata.long, x = "scaled",
  #             add = "mean", rug = TRUE,
  #             color = "Treatment", fill = "Treatment",
  #             facet.by = 'transcript') +
  #   scale_fill_manual(name = "", values = rev(p_colors)) +
  #   scale_color_manual(name = "", values = rev(p_colors)) +
  #   ggtitle("scaled") +
  #   facet_wrap(~transcript, ncol = 4, scales = "free")
  # print(p2)
  
  
  cat(crayon::red('####  ####  \nDistribution tests\n####  ####  \n'))
  
    distr_tests = mydata.long %>%
    dplyr::group_by(transcript) %>%
    dplyr::group_map(~ {
      x = na.omit(.x$measurement)
      if (length(x) < 3) {
        return(tibble::tibble(
          Shapiro_Wilk = NA_real_,
          Anderson_Darling = NA_real_,
          Lilliefors_KS = NA_real_,
          Jarque_Bera = NA_real_,
          DAgostino_Skewness = NA_real_,
          n = length(x),
          transcript = .y$transcript
        ))
      }
      tibble::tibble(
        Shapiro_Wilk = tryCatch(shapiro.test(x)$p.value, error = function(e) NA_real_),
        Anderson_Darling = tryCatch(nortest::ad.test(x)$p.value, error = function(e) NA_real_),
        Lilliefors_KS = tryCatch(nortest::lillie.test(x)$p.value, error = function(e) NA_real_),
        Jarque_Bera = tryCatch(tseries::jarque.bera.test(x)$p.value, error = function(e) NA_real_),
        DAgostino_Skewness = tryCatch(moments::agostino.test(x)$p.value, error = function(e) NA_real_),
        n = length(x),
        transcript = .y$transcript
      )
    }) %>%
    dplyr::bind_rows() %>%
    dplyr::ungroup() %>%
    dplyr::mutate(
      Shapiro_Wilk_BH = stats::p.adjust(Shapiro_Wilk, method = "BH"),
      Anderson_Darling_BH = stats::p.adjust(Anderson_Darling, method = "BH"),
      Lilliefors_KS_BH = stats::p.adjust(Lilliefors_KS, method = "BH"),
      Jarque_Bera_BH = stats::p.adjust(Jarque_Bera, method = "BH"),
      DAgostino_Skewness_BH = stats::p.adjust(DAgostino_Skewness, method = "BH")
    )

  

  combined_results = distr_tests[grep("^n$|transcript|_BH", colnames(distr_tests))]
  print(combined_results)

  
  

  cat(crayon::red("####  ####  \nQuantile-Quantile plots\n####  ####  \n"))
  qq_plot_list = mydata.long %>%
    dplyr::group_by(transcript) %>%
    dplyr::group_map(~ {
      ggpubr::ggqqplot(.x$measurement) +
        ggplot2::ggtitle(paste("Raw:", .y$transcript))
    })
  
  resid_qq_list = mydata.long %>%
    dplyr::group_by(transcript) %>%
    dplyr::group_map(~ {
      model = lm(measurement ~ Treatment, data = .x)
      resids = resid(model)
      ggpubr::ggqqplot(resids) +
        ggplot2::ggtitle(paste("Residuals:", .y$transcript))
    })
  

  resid_density_list = mydata.long %>%
    dplyr::group_by(transcript) %>%
    dplyr::group_map(~ {
      model = lm(measurement ~ Treatment, data = .x)
      resids = scale(resid(model))  # standardize residuals
      df = data.frame(resids = as.numeric(resids))
      ggplot(df, aes(x = resids)) +
        geom_density(fill = "steelblue", alpha = 0.6) +
        stat_function(fun = dnorm, color = "red", linetype = "dashed") +
        ggtitle(paste("Residuals:", .y$transcript)) +
        theme_minimal()
    })
  
  # Arrange and print all plots
  qq_arranged = ggpubr::ggarrange(plotlist = qq_plot_list, ncol = 4, nrow = ceiling(length(qq_plot_list) / 4))
  resid_qq_arranged = ggpubr::ggarrange(plotlist = resid_qq_list, ncol = 4, nrow = ceiling(length(resid_qq_list) / 4))
  resid_hist_arranged = ggpubr::ggarrange(plotlist = resid_density_list, ncol = 4, nrow = ceiling(length(resid_density_list) / 4))
  
  print(qq_arranged)
  print(resid_qq_arranged)
  print(resid_hist_arranged)
  
  # Return invisibly
  invisible(list(qq = qq_arranged, resid_qq = resid_qq_arranged, resid_hist = resid_hist_arranged))



  
  cat(crayon::red("####  ####  \nTest for homogeneity of variance across groups\n####  ####  \n"))
  
  cat(crayon::red("####  ####  \nLevene\n####  ####  \n"))
  # Levene's test on measurement
  lev_meas = mydata.long %>%
    dplyr::group_by(transcript) %>%
    rstatix::levene_test(measurement ~ Treatment, center = "mean")
  print(lev_meas)
  
  cat(crayon::red("####  ####  \nBrown-Forsythe\n####  ####  \n"))
  # center = "median" switches Levene’s test to the Brown-Forsythe variant
  bf_meas = mydata.long %>%
    dplyr::group_by(transcript) %>%
    rstatix::levene_test(measurement ~ Treatment, center = "median")
  print(bf_meas)
  
  cat(crayon::red("####  ####  \nFligner\n####  ####  \n"))
  fk_meas = mydata.long %>%
    dplyr::group_by(transcript) %>%
    dplyr::group_modify(~ broom::tidy(stats::fligner.test(measurement ~ Treatment, data = .x)))
  print(fk_meas)


  

  cat(crayon::red("####  ####  \nWilcoxon effect size\n####  ####  \n"))
  
  # Wilcoxon effect size on measurement
  eff_meas = mydata.long %>%
    dplyr::group_by(transcript) %>%
    rstatix::wilcox_effsize(measurement ~ Treatment)
  print(eff_meas)
  
  # Wilcoxon effect size on scaled
  # eff_scaled = mydata.long %>%
  #   dplyr::group_by(transcript) %>%
  #   rstatix::wilcox_effsize(scaled ~ Treatment)
  # print(eff_scaled)
  
  cat(crayon::red("####  ####  \nCohen's d Measure of Effect Size\n####  ####  \n"))
  
  # Cohen's d on measurement
  coh_meas = mydata.long %>%
    dplyr::group_by(transcript) %>%
    rstatix::cohens_d(measurement ~ Treatment, paired = FALSE)
  print(coh_meas)
  

  invisible(list(
    p_measurement = p1,
    distr_measurement = combined_results,
    levene_measurement = lev_meas,
    brownforsythe_measurement = bf_meas,
    fligner_measurement = fk_meas,
    wilcoxon_measurement = eff_meas,
    cohensd_measurement = coh_meas
  ))
}

3.6 ggplot limits

make_scale_params <- function(maxval) {
  upper = ifelse(maxval > 2, maxval + 2, maxval + 0.1) # 0.5
  step = dplyr::case_when(
    maxval >= 50 ~ 25,
    maxval >= 30 ~ 10,
    maxval >= 15 ~ 5,
    maxval >= 10 ~ 2,
    maxval >  2  ~ 1,
    TRUE         ~ 0.5
  )
  upper_round = ceiling(upper / step) * step
  breaks = seq(0, upper_round, by = step)
  list(limits = c(0, upper_round), breaks = breaks)
}

build_y_scales_for <- function(names_vec, max_per_transcript) {
  max_per_transcript %>%
    dplyr::filter(transcript %in% names_vec) %>%
    purrr::pmap(function(transcript, max_scaled) {
      params = make_scale_params(max_scaled)
      lhs = rlang::expr(transcript == !!transcript)
      rhs = rlang::expr(ggplot2::scale_y_continuous(limits = !!params$limits, breaks = !!params$breaks))
      rlang::new_formula(lhs, rhs)
    })
}

3.7 test and plot

test_and_plot <- function(data_long_raw
                          , myorder
                          , pal
                          , what
                          , plot_gene_expression_func
                          , groupvars = c("Treatment", "transcript")
                          , y_scales1
                          , y_scales2
                          ) {

  
  mydata.long = dplyr::as_tibble(data.table::data.table(data_long_raw))
  mydata.long$transcript = factor(mydata.long$transcript, levels = myorder)
  mydata.long = mydata.long %>% dplyr::arrange(factor(transcript, levels = myorder))
  
  mydata.long.SE = summary_stats(mydata.long, measurevar = "scaled", groupvars = groupvars)
  
  results = dens_and_effect(mydata.long, p_colors = pal)
  
  cat(what, file = fr, append = TRUE, sep = "\n")
  
  cat("\nDistribution tests", file = fr, append = TRUE, sep = "\n")
  header = paste(colnames(results$distr_measurement), collapse = "\t")
  cat(header, file = fr, append = TRUE, sep = "\n")
  output_text = apply(results$distr_measurement, 1, function(row) paste(row, collapse = "\t"))
  cat(output_text, file = fr, append = TRUE, sep = "\n")
  
  cat("\nLevene's test for homogeneity of variance across groups", file = fr, append = TRUE, sep = "\n")
  header = paste(colnames(results$levene_measurement), collapse = "\t")
  cat(header, file = fr, append = TRUE, sep = "\n")
  output_text = apply(results$levene_measurement, 1, function(row) paste(row, collapse = "\t"))
  cat(output_text, file = fr, append = TRUE, sep = "\n")
  cat("\nBrown-Forsythe (Levene with median)", file = fr, append = TRUE, sep = "\n")
  header = paste(colnames(results$brownforsythe_measurement), collapse = "\t")
  cat(header, file = fr, append = TRUE, sep = "\n")
  output_text = apply(results$brownforsythe_measurement, 1, function(row) paste(row, collapse = "\t"))
  cat(output_text, file = fr, append = TRUE, sep = "\n")
  cat("\nFligner-Killeen test of homogeneity of variances", file = fr, append = TRUE, sep = "\n")
  header = paste(colnames(results$fligner_measurement), collapse = "\t")
  cat(header, file = fr, append = TRUE, sep = "\n")
  output_text = apply(results$fligner_measurement, 1, function(row) paste(row, collapse = "\t"))
  cat(output_text, file = fr, append = TRUE, sep = "\n")
  
  cat("\nWilcoxon effect size", file = fr, append = TRUE, sep = "\n")
  header = paste(colnames(results$wilcoxon_measurement), collapse = "\t")
  cat(header, file = fr, append = TRUE, sep = "\n")
  output_text = apply(results$wilcoxon_measurement, 1, function(row) paste(row, collapse = "\t"))
  cat(output_text, file = fr, append = TRUE, sep = "\n")
  cat("\nCohen's d Measure of Effect Size", file = fr, append = TRUE, sep = "\n")
  header = paste(colnames(results$cohensd_measurement), collapse = "\t")
  cat(header, file = fr, append = TRUE, sep = "\n")
  output_text = apply(results$cohensd_measurement, 1, function(row) paste(row, collapse = "\t"))
  cat(output_text, file = fr, append = TRUE, sep = "\n")
  
  cat("", file = fr, append = TRUE, sep = "\n")


  
  temp = perm_test_by_transcript(mydata.long, measurevar = "measurement", groupvar = "Treatment")
  temp$transcript = rownames(temp)
  perm = tidyr::gather(temp, contrast, perm.p, colnames(temp)[1]:colnames(temp)[ncol(temp)-1], factor_key=TRUE)
  perm$group1 = gsub(' vs.*', '', perm$contrast)
  perm$group2 = gsub('.* vs ', '', perm$contrast)
  perm$perm.p.adj = p.adjust(perm$perm.p, method = 'BH')
  perm$perm.p.adj.signif = 'ns'
  perm$perm.p.adj.signif[perm$perm.p.adj < 0.0001] = '**'
  perm$perm.p.adj.signif[perm$perm.p.adj < 0.001] = '**'
  perm$perm.p.adj.signif[perm$perm.p.adj < 0.05] = '*'
  
  
  stat.test = perm %>%
    dplyr::select(transcript, group1, group2, perm.p, perm.p.adj, perm.p.adj.signif) %>%
    dplyr::rename(
      p = perm.p,
      p.adj = perm.p.adj,
      p.adj.signif = perm.p.adj.signif
    )
  
  # Compute y-position using both mean ± se and scaled values
  y_pos_df = mydata.long.SE %>%
    dplyr::group_by(transcript) %>%
    dplyr::summarise(
      max_mean_se = max(mean + se, na.rm = TRUE),
      min_mean_se = min(mean - se, na.rm = TRUE),
      .groups = "drop"
    ) %>%
    dplyr::left_join(
      mydata.long %>%
        dplyr::group_by(transcript) %>%
        dplyr::summarise(
          max_scaled = max(scaled, na.rm = TRUE),
          min_scaled = min(scaled, na.rm = TRUE),
          .groups = "drop"
        ),
      by = "transcript"
    ) %>%
    dplyr::mutate(
      plot_max = pmax(max_mean_se, max_scaled, na.rm = TRUE),
      plot_min = pmin(min_mean_se, min_scaled, na.rm = TRUE),
      plot_range = pmax(plot_max - plot_min, 1e-6)
    ) %>%
    dplyr::select(transcript, plot_max, plot_range)
  
  # Map group names to numeric x positions
  pos_map = mydata.long %>%
    dplyr::distinct(transcript, Treatment) %>%
    dplyr::mutate(xpos = as.numeric(factor(Treatment, levels = levels(mydata.long$Treatment))))
  
  # Attach y-position and x-position info to stat.test
  stat.test = stat.test %>%
    dplyr::left_join(y_pos_df, by = "transcript") %>%
    dplyr::left_join(pos_map %>% dplyr::rename(group1 = Treatment, xmin = xpos), by = c("transcript", "group1")) %>%
    dplyr::left_join(pos_map %>% dplyr::rename(group2 = Treatment, xmax = xpos), by = c("transcript", "group2")) %>%
    dplyr::mutate(
      xmin = dplyr::if_else(is.na(xmin), as.numeric(factor(group1, levels = levels(mydata.long$Treatment))), xmin),
      xmax = dplyr::if_else(is.na(xmax), as.numeric(factor(group2, levels = levels(mydata.long$Treatment))), xmax)
    )
  
  # Compute distinct y.position per comparison for significant results
  stat.test.sig = stat.test %>%
    dplyr::filter(!is.na(p.adj) & p.adj <= 0.05) %>%
    dplyr::group_by(transcript) %>%
    dplyr::arrange(xmin, xmax, .by_group = TRUE) %>%
    dplyr::mutate(
      base_y = plot_max + 0.03 * plot_range,
      inc = 0.05 * plot_range,
      y.position = base_y + (dplyr::row_number() - 1) * inc
    ) %>%
    dplyr::ungroup() %>%
    dplyr::select(transcript, group1, group2, p.adj, p.adj.signif, y.position, xmin, xmax)

  
  group1 =  sapply(y_scales1, function(x) {
    s = deparse(x)
    s_full = paste(s, collapse = " ")
    sub('.*transcript == *"([^"]+)".*', '\\1', s_full)
  })
  group2 =  sapply(y_scales2, function(x) {
    s = deparse(x)
    s_full = paste(s, collapse = " ")
    sub('.*transcript == *"([^"]+)".*', '\\1', s_full)
  })


  p1 = plot_gene_expression_func(
    data.SE = mydata.long.SE,
    data.long = mydata.long,
    stat.test.sig = stat.test.sig,
    transcripts_excl = group2,
    facet_cols = 2,
    color_values = pal,
    plot_title = what,
    y_scales = y_scales1,
    dodge_width = 0.8
  )
  print(p1)
  
  p2 = plot_gene_expression_func(
    data.SE = mydata.long.SE,
    data.long = mydata.long,
    stat.test.sig = stat.test.sig,
    transcripts_excl = group1,
    facet_cols = 6,
    color_values = pal,
    plot_title = what,
    y_scales = y_scales2,
    dodge_width = 0.8
  )
  print(p2)
  
  return(list(plot1 = p1, plot2 = p2, stat.test = stat.test))
}

4 Desiree

4.1 input

fp = file.path('..', "input")
list.files(fp)
## [1] "~$Table_combinedInputs.xlsx" "README.md"                  
## [3] "Table_combinedInputs.xlsx"
fn = 'Table_combinedInputs.xlsx'
wb = openxlsx::loadWorkbook(file.path(fp, fn))
sheet_names = names(wb)
sheet_name = grep("Desiree_Rywal", sheet_names, value = TRUE)


myTable = openxlsx::read.xlsx(xlsxFile = file.path(fp, fn),
                            sheet = sheet_name,
                            startRow = 10,
                            colNames = TRUE,
                            rowNames = FALSE,
                            detectDates = FALSE,
                            skipEmptyRows = TRUE,
                            skipEmptyCols = TRUE,
                            rows = NULL,
                            cols = NULL,
                            check.names = FALSE,
                            sep.names = ".",
                            namedRegion = NULL,
                            na.strings = "NA",
                            fillMergedCells = FALSE)

myTable$Genotype = trimws(myTable$Genotype)
Desiree = myTable[myTable$Genotype == 'Desiree', ]

data.table::setDT(Desiree)
Desiree[, SampleID := NULL]
Desiree[, Genotype := NULL]

Desiree$Tissue = as.factor(trimws(Desiree$Tissue))
Desiree$Treatment = factor(trimws(Desiree$Treatment), levels = c("noninoculated", "inoculated"))

4.2 params

my_dir = file.path("..", "reports", "Desiree_Rywal")
if (!dir.exists(my_dir)) {
  dir.create(my_dir)
}
fr = file.path('..', 'output',  'Desiree_Rywal_stat.txt')
file.create(fr)
## [1] TRUE
groupvars = c("Tissue", "Treatment", "transcript")
measurevar = 'measurement'
dodge = position_dodge(width = 0.5)

myorder = c("StRBOHD", "StPR1B", "StCPI8", "StCAB", "StBGLU2", "StHSP70", "StPti5", "St13-LOX")
group1_names = c("StPti5", "St13-LOX")

4.3 results

cat("Desiree", file = fr, append = TRUE, sep = "\n")
cat("", file = fr, append = TRUE, sep = "\n")


run_analysis_for_strain <- function(data
                                  , strain
                                  , myorder
                                  , pal
                                  , groupvars
                                  , measurevar
                                  , my_dir
                                  , plot_gene_expression
                                  , test_and_plot) {


  strains = trimws(unique(na.omit(data$Strain)))
  temp = data[grep(setdiff(strains, strain), data$Strain, invert = TRUE), ]
  temp[, Strain := NULL]


  shoots = process_data(temp, groupvars, measurevar, scale_treatment = "noninoculated")$shoots
  roots = process_data(temp, groupvars, measurevar, scale_treatment = "noninoculated")$roots
  
  max_per_transcript = shoots %>%
    dplyr::group_by(transcript) %>%
    dplyr::summarise(max_scaled = max(scaled, na.rm = TRUE), .groups = "drop")
  group2_names = setdiff(max_per_transcript$transcript, group1_names)
  y_scales1 = build_y_scales_for(group1_names, max_per_transcript)
  y_scales2 = build_y_scales_for(group2_names, max_per_transcript)
  max_per_transcript = roots %>%
    dplyr::group_by(transcript) %>%
    dplyr::summarise(max_scaled = max(scaled, na.rm = TRUE), .groups = "drop")
  y_scales3 = build_y_scales_for(group1_names, max_per_transcript)
  y_scales4 = build_y_scales_for(group2_names, max_per_transcript)


  

  result_shoots = test_and_plot(data_long_raw = shoots,
                                myorder = myorder,
                                pal = pal,
                                what = paste("shoots", strain),
                                plot_gene_expression_func = plot_gene_expression,
                                groupvars = groupvars,
                                y_scales1 = y_scales1,
                                y_scales2 = y_scales2)
  res = result_shoots$stat.test[, grep("transcript|group1|group2|^p$|p\\.", colnames(result_shoots$stat.test))]
  print(res[res$p.adj.signif != 'ns', ])
  
  cat("", file = fr, append = TRUE, sep = "\n")
  output_text = "permutational t-test"
  cat(output_text, file = fr, append = TRUE, sep = "\n")
  header = paste(colnames(res), collapse = "\t")
  cat(header, file = fr, append = TRUE, sep = "\n")
  output_text = apply(as.data.frame(tibble::as_tibble(res)), 1, function(row) paste(row, collapse = "\t"))
  cat(output_text, file = fr, append = TRUE, sep = "\n")
  cat("", file = fr, append = TRUE, sep = "\n")
  cat("", file = fr, append = TRUE, sep = "\n")
  
  
  ggsave(filename = file.path(my_dir, paste0("Desiree_shoots.", gsub(" ", "", strain), "_1.pdf")),
         plot = result_shoots$plot1, device = pdf, width = 3, height = 8, units = "in", dpi = 900)
  ggsave(filename = file.path(my_dir, paste0("Desiree_shoots.", gsub(" ", "", strain), "_1.svg")),
         plot = result_shoots$plot1, device = svg, width = 3, height = 8, units = "in", dpi = 900)
    
  ggsave(filename = file.path(my_dir, paste0("Desiree_shoots.", gsub(" ", "", strain), "_2.pdf")),
         plot = result_shoots$plot2, device = pdf, width = 9, height = 8, units = "in", dpi = 900)
  ggsave(filename = file.path(my_dir, paste0("Desiree_shoots.", gsub(" ", "", strain), "_2.svg")),
         plot = result_shoots$plot2, device = svg, width = 9, height = 8, units = "in", dpi = 900)
  

  result_roots = test_and_plot(data_long_raw = roots,
                               myorder = myorder,
                               pal = pal,
                               what = paste("roots", strain),
                               plot_gene_expression_func = plot_gene_expression,
                               groupvars = groupvars,
                               y_scales1 = y_scales3,
                               y_scales2 = y_scales4)
  res = result_roots$stat.test[, grep("transcript|group1|group2|^p$|p\\.", colnames(result_roots$stat.test))]
  print(res[res$p.adj.signif != 'ns', ])
  
  cat("", file = fr, append = TRUE, sep = "\n")
  output_text = "permutational t-test"
  cat(output_text, file = fr, append = TRUE, sep = "\n")
  header = paste(colnames(res), collapse = "\t")
  cat(header, file = fr, append = TRUE, sep = "\n")
  output_text = apply(as.data.frame(tibble::as_tibble(res)), 1, function(row) paste(row, collapse = "\t"))
  cat(output_text, file = fr, append = TRUE, sep = "\n")
  cat("", file = fr, append = TRUE, sep = "\n")
  cat("", file = fr, append = TRUE, sep = "\n")
  
  
  ggsave(filename = file.path(my_dir, paste0("Desiree_roots.", gsub(" ", ".", strain), "_1.pdf")),
         plot = result_roots$plot1, device = pdf, width = 3, height = 8, units = "in", dpi = 900)
  ggsave(filename = file.path(my_dir, paste0("Desiree_roots.", gsub(" ", ".", strain), "_1.svg")),
         plot = result_roots$plot1, device = svg, width = 3, height = 8, units = "in", dpi = 900)
  
  ggsave(filename = file.path(my_dir, paste0("Desiree_roots.", gsub(" ", ".", strain), "_2.pdf")),
         plot = result_roots$plot2, device = pdf, width = 9, height = 8, units = "in", dpi = 900)
  ggsave(filename = file.path(my_dir, paste0("Desiree_roots.", gsub(" ", ".", strain), "_2.svg")),
         plot = result_roots$plot2, device = svg, width = 9, height = 8, units = "in", dpi = 900)
  

  list(shoots = result_shoots, roots = result_roots)
}


results_PS216 = run_analysis_for_strain(data = Desiree
                                   , strain = "PS-216"
                                   , myorder
                                   , pal
                                   , groupvars
                                   , measurevar
                                   , my_dir
                                   , plot_gene_expression
                                   , test_and_plot)

## ####  ####  
## Distribution tests
## ####  ####
## # A tibble: 8 × 7
##       n transcript Shapiro_Wilk_BH Anderson_Darling_BH Lilliefors_KS_BH
##   <int> <fct>                <dbl>               <dbl>            <dbl>
## 1    12 StRBOHD           0.0222              0.0554            0.471  
## 2    12 StPR1B            0.495               0.585             0.604  
## 3    10 StCPI8            0.132               0.139             0.239  
## 4     9 StCAB             0.115               0.112             0.246  
## 5    11 StBGLU2           0.000110            0.000245          0.00649
## 6    12 StHSP70           0.495               0.585             0.578  
## 7    12 StPti5            0.0158              0.0141            0.0298 
## 8    12 St13-LOX          0.853               0.796             0.604  
## # ℹ 2 more variables: Jarque_Bera_BH <dbl>, DAgostino_Skewness_BH <dbl>
## ####  ####  
## Quantile-Quantile plots
## ####  ####

## ####  ####  
## Test for homogeneity of variance across groups
## ####  ####  
## ####  ####  
## Levene
## ####  ####  
## # A tibble: 8 × 5
##   transcript   df1   df2 statistic       p
##   <fct>      <int> <int>     <dbl>   <dbl>
## 1 StRBOHD        1    10     3.44  0.0935 
## 2 StPR1B         1    10     0.925 0.359  
## 3 StCPI8         1     8     5.29  0.0505 
## 4 StCAB          1     7     1.97  0.203  
## 5 StBGLU2        1     9     2.66  0.137  
## 6 StHSP70        1    10    10.1   0.00997
## 7 StPti5         1    10     9.19  0.0127 
## 8 St13-LOX       1    10     1.98  0.190  
## ####  ####  
## Brown-Forsythe
## ####  ####  
## # A tibble: 8 × 5
##   transcript   df1   df2 statistic      p
##   <fct>      <int> <int>     <dbl>  <dbl>
## 1 StRBOHD        1    10     2.15  0.173 
## 2 StPR1B         1    10     0.769 0.401 
## 3 StCPI8         1     8     4.30  0.0719
## 4 StCAB          1     7     1.34  0.286 
## 5 StBGLU2        1     9     0.639 0.444 
## 6 StHSP70        1    10     1.02  0.335 
## 7 StPti5         1    10     3.33  0.0981
## 8 St13-LOX       1    10     0.716 0.417 
## ####  ####  
## Fligner
## ####  ####  
## # A tibble: 8 × 5
## # Groups:   transcript [8]
##   transcript statistic p.value parameter method                                 
##   <fct>          <dbl>   <dbl>     <dbl> <chr>                                  
## 1 StRBOHD      4.16     0.0414         1 Fligner-Killeen test of homogeneity of…
## 2 StPR1B       0.587    0.444          1 Fligner-Killeen test of homogeneity of…
## 3 StCPI8       2.52     0.113          1 Fligner-Killeen test of homogeneity of…
## 4 StCAB        1.31     0.252          1 Fligner-Killeen test of homogeneity of…
## 5 StBGLU2      0.0112   0.916          1 Fligner-Killeen test of homogeneity of…
## 6 StHSP70      0.00257  0.960          1 Fligner-Killeen test of homogeneity of…
## 7 StPti5       5.29     0.0215         1 Fligner-Killeen test of homogeneity of…
## 8 St13-LOX     0.254    0.614          1 Fligner-Killeen test of homogeneity of…
## ####  ####  
## Wilcoxon effect size
## ####  ####  
## # A tibble: 8 × 8
##   .y.         group1        group2     effsize transcript    n1    n2 magnitude
## * <chr>       <chr>         <chr>        <dbl> <fct>      <int> <int> <ord>    
## 1 measurement noninoculated inoculated   0.832 StRBOHD        6     6 large    
## 2 measurement noninoculated inoculated   0.139 StPR1B         6     6 small    
## 3 measurement noninoculated inoculated   0.826 StCPI8         5     5 large    
## 4 measurement noninoculated inoculated   0.816 StCAB          5     4 large    
## 5 measurement noninoculated inoculated   0.220 StBGLU2        5     6 small    
## 6 measurement noninoculated inoculated   0.277 StHSP70        6     6 small    
## 7 measurement noninoculated inoculated   0.786 StPti5         6     6 large    
## 8 measurement noninoculated inoculated   0.832 St13-LOX       6     6 large    
## ####  ####  
## Cohen's d Measure of Effect Size
## ####  ####  
## # A tibble: 8 × 8
##   .y.         group1        group2     effsize transcript    n1    n2 magnitude
## * <chr>       <chr>         <chr>        <dbl> <fct>      <int> <int> <ord>    
## 1 measurement noninoculated inoculated   1.81  StRBOHD        6     6 large    
## 2 measurement noninoculated inoculated  -0.386 StPR1B         6     6 small    
## 3 measurement noninoculated inoculated  -4.61  StCPI8         5     5 large    
## 4 measurement noninoculated inoculated   7.35  StCAB          5     4 large    
## 5 measurement noninoculated inoculated  -0.595 StBGLU2        5     6 moderate 
## 6 measurement noninoculated inoculated  -0.215 StHSP70        6     6 small    
## 7 measurement noninoculated inoculated  -1.49  StPti5         6     6 large    
## 8 measurement noninoculated inoculated  -2.82  St13-LOX       6     6 large

##   transcript        group1     group2            p        p.adj p.adj.signif
## 1    StRBOHD noninoculated inoculated 0.0002442599 0.0004885198            *
## 3     StCPI8 noninoculated inoculated 0.0002442599 0.0004885198            *
## 4      StCAB noninoculated inoculated 0.0002442599 0.0004885198            *
## 7     StPti5 noninoculated inoculated 0.0036638984 0.0058622374            *
## 8   St13-LOX noninoculated inoculated 0.0002442599 0.0004885198            *

## ####  ####  
## Distribution tests
## ####  ####  
## # A tibble: 8 × 7
##       n transcript Shapiro_Wilk_BH Anderson_Darling_BH Lilliefors_KS_BH
##   <int> <fct>                <dbl>               <dbl>            <dbl>
## 1    12 StRBOHD           0.0785            0.0784           0.0433    
## 2    12 StPR1B            0.000103          0.00000113       0.00000420
## 3    12 StCPI8            0.558             0.670            0.769     
## 4     9 StCAB             0.558             0.670            0.769     
## 5    12 StBGLU2           0.0785            0.0830           0.0433    
## 6    12 StHSP70           0.809             0.827            0.769     
## 7    11 StPti5            0.0102            0.0139           0.00702   
## 8    12 St13-LOX          0.657             0.768            0.769     
## # ℹ 2 more variables: Jarque_Bera_BH <dbl>, DAgostino_Skewness_BH <dbl>
## ####  ####  
## Quantile-Quantile plots
## ####  ####

## ####  ####  
## Test for homogeneity of variance across groups
## ####  ####  
## ####  ####  
## Levene
## ####  ####  
## # A tibble: 8 × 5
##   transcript   df1   df2 statistic      p
##   <fct>      <int> <int>     <dbl>  <dbl>
## 1 StRBOHD        1    10     5.58  0.0399
## 2 StPR1B         1    10     4.63  0.0569
## 3 StCPI8         1    10     0.246 0.630 
## 4 StCAB          1     7     2.28  0.175 
## 5 StBGLU2        1    10     4.77  0.0539
## 6 StHSP70        1    10     1.58  0.238 
## 7 StPti5         1     9     7.91  0.0203
## 8 St13-LOX       1    10     1.59  0.236 
## ####  ####  
## Brown-Forsythe
## ####  ####  
## # A tibble: 8 × 5
##   transcript   df1   df2 statistic      p
##   <fct>      <int> <int>     <dbl>  <dbl>
## 1 StRBOHD        1    10     4.02  0.0729
## 2 StPR1B         1    10     1.38  0.268 
## 3 StCPI8         1    10     0.354 0.565 
## 4 StCAB          1     7     2.48  0.159 
## 5 StBGLU2        1    10     4.18  0.0682
## 6 StHSP70        1    10     1.30  0.281 
## 7 StPti5         1     9     6.44  0.0318
## 8 St13-LOX       1    10     1.45  0.257 
## ####  ####  
## Fligner
## ####  ####  
## # A tibble: 8 × 5
## # Groups:   transcript [8]
##   transcript statistic p.value parameter method                                 
##   <fct>          <dbl>   <dbl>     <dbl> <chr>                                  
## 1 StRBOHD        2.99  0.0840          1 Fligner-Killeen test of homogeneity of…
## 2 StPR1B         1.14  0.285           1 Fligner-Killeen test of homogeneity of…
## 3 StCPI8         0.521 0.470           1 Fligner-Killeen test of homogeneity of…
## 4 StCAB          2.89  0.0890          1 Fligner-Killeen test of homogeneity of…
## 5 StBGLU2        3.41  0.0649          1 Fligner-Killeen test of homogeneity of…
## 6 StHSP70        2.27  0.132           1 Fligner-Killeen test of homogeneity of…
## 7 StPti5         6.64  0.00998         1 Fligner-Killeen test of homogeneity of…
## 8 St13-LOX       0.784 0.376           1 Fligner-Killeen test of homogeneity of…
## ####  ####  
## Wilcoxon effect size
## ####  ####  
## # A tibble: 8 × 8
##   .y.         group1        group2     effsize transcript    n1    n2 magnitude
## * <chr>       <chr>         <chr>        <dbl> <fct>      <int> <int> <ord>    
## 1 measurement noninoculated inoculated   0.740 StRBOHD        6     6 large    
## 2 measurement noninoculated inoculated   0.233 StPR1B         6     6 small    
## 3 measurement noninoculated inoculated   0     StCPI8         6     6 small    
## 4 measurement noninoculated inoculated   0.602 StCAB          6     3 large    
## 5 measurement noninoculated inoculated   0.370 StBGLU2        6     6 moderate 
## 6 measurement noninoculated inoculated   0.324 StHSP70        6     6 moderate 
## 7 measurement noninoculated inoculated   0.330 StPti5         5     6 moderate 
## 8 measurement noninoculated inoculated   0.832 St13-LOX       6     6 large    
## ####  ####  
## Cohen's d Measure of Effect Size
## ####  ####  
## # A tibble: 8 × 8
##   .y.         group1        group2     effsize transcript    n1    n2 magnitude 
## * <chr>       <chr>         <chr>        <dbl> <fct>      <int> <int> <ord>     
## 1 measurement noninoculated inoculated  2.37   StRBOHD        6     6 large     
## 2 measurement noninoculated inoculated  0.647  StPR1B         6     6 moderate  
## 3 measurement noninoculated inoculated -0.0549 StCPI8         6     6 negligible
## 4 measurement noninoculated inoculated  1.56   StCAB          6     3 large     
## 5 measurement noninoculated inoculated  0.850  StBGLU2        6     6 large     
## 6 measurement noninoculated inoculated  0.669  StHSP70        6     6 moderate  
## 7 measurement noninoculated inoculated -1.00   StPti5         5     6 large     
## 8 measurement noninoculated inoculated -2.64   St13-LOX       6     6 large

##   transcript        group1     group2            p       p.adj p.adj.signif
## 1    StRBOHD noninoculated inoculated 0.0031753786 0.012701514            *
## 8   St13-LOX noninoculated inoculated 0.0002442599 0.001954079            *

results_PS218 = run_analysis_for_strain(data = Desiree
                                   , strain = "PS-218"
                                   , myorder
                                   , pal
                                   , groupvars
                                   , measurevar
                                   , my_dir
                                   , plot_gene_expression
                                   , test_and_plot)

## ####  ####  
## Distribution tests
## ####  ####  
## # A tibble: 8 × 7
##       n transcript Shapiro_Wilk_BH Anderson_Darling_BH Lilliefors_KS_BH
##   <int> <fct>                <dbl>               <dbl>            <dbl>
## 1    12 StRBOHD            0.0286              0.0580            0.436 
## 2    12 StPR1B             0.0253              0.0358            0.117 
## 3    11 StCPI8             0.00614             0.00812           0.0220
## 4    11 StCAB              0.0452              0.0443            0.0458
## 5    11 StBGLU2            0.100               0.101             0.117 
## 6    12 StHSP70            0.0253              0.0120            0.0220
## 7    11 StPti5             0.0103              0.00812           0.0220
## 8    12 St13-LOX           0.289               0.209             0.117 
## # ℹ 2 more variables: Jarque_Bera_BH <dbl>, DAgostino_Skewness_BH <dbl>
## ####  ####  
## Quantile-Quantile plots
## ####  ####

## ####  ####  
## Test for homogeneity of variance across groups
## ####  ####  
## ####  ####  
## Levene
## ####  ####  
## # A tibble: 8 × 5
##   transcript   df1   df2 statistic       p
##   <fct>      <int> <int>     <dbl>   <dbl>
## 1 StRBOHD        1    10     2.07  0.181  
## 2 StPR1B         1    10     0.843 0.380  
## 3 StCPI8         1     9     8.34  0.0179 
## 4 StCAB          1     9     2.03  0.188  
## 5 StBGLU2        1     9     0.484 0.504  
## 6 StHSP70        1    10     1.06  0.328  
## 7 StPti5         1     9    17.9   0.00219
## 8 St13-LOX       1    10     6.47  0.0292 
## ####  ####  
## Brown-Forsythe
## ####  ####  
## # A tibble: 8 × 5
##   transcript   df1   df2 statistic       p
##   <fct>      <int> <int>     <dbl>   <dbl>
## 1 StRBOHD        1    10    1.32   0.277  
## 2 StPR1B         1    10    0.0573 0.816  
## 3 StCPI8         1     9    2.22   0.171  
## 4 StCAB          1     9    1.43   0.263  
## 5 StBGLU2        1     9    0.0630 0.808  
## 6 StHSP70        1    10    0.113  0.743  
## 7 StPti5         1     9   11.6    0.00771
## 8 St13-LOX       1    10    5.91   0.0354 
## ####  ####  
## Fligner
## ####  ####  
## # A tibble: 8 × 5
## # Groups:   transcript [8]
##   transcript statistic p.value parameter method                                 
##   <fct>          <dbl>   <dbl>     <dbl> <chr>                                  
## 1 StRBOHD      1.41     0.234          1 Fligner-Killeen test of homogeneity of…
## 2 StPR1B       0.327    0.567          1 Fligner-Killeen test of homogeneity of…
## 3 StCPI8       2.89     0.0890         1 Fligner-Killeen test of homogeneity of…
## 4 StCAB        0.640    0.424          1 Fligner-Killeen test of homogeneity of…
## 5 StBGLU2      0.0111   0.916          1 Fligner-Killeen test of homogeneity of…
## 6 StHSP70      0.00469  0.945          1 Fligner-Killeen test of homogeneity of…
## 7 StPti5       3.36     0.0666         1 Fligner-Killeen test of homogeneity of…
## 8 St13-LOX     5.30     0.0213         1 Fligner-Killeen test of homogeneity of…
## ####  ####  
## Wilcoxon effect size
## ####  ####  
## # A tibble: 8 × 8
##   .y.         group1        group2     effsize transcript    n1    n2 magnitude
## * <chr>       <chr>         <chr>        <dbl> <fct>      <int> <int> <ord>    
## 1 measurement noninoculated inoculated   0.740 StRBOHD        6     6 large    
## 2 measurement noninoculated inoculated   0.139 StPR1B         6     6 small    
## 3 measurement noninoculated inoculated   0.826 StCPI8         5     6 large    
## 4 measurement noninoculated inoculated   0.826 StCAB          5     6 large    
## 5 measurement noninoculated inoculated   0.165 StBGLU2        5     6 small    
## 6 measurement noninoculated inoculated   0.231 StHSP70        6     6 small    
## 7 measurement noninoculated inoculated   0.826 StPti5         6     5 large    
## 8 measurement noninoculated inoculated   0.832 St13-LOX       6     6 large    
## ####  ####  
## Cohen's d Measure of Effect Size
## ####  ####  
## # A tibble: 8 × 8
##   .y.         group1        group2     effsize transcript    n1    n2 magnitude 
## * <chr>       <chr>         <chr>        <dbl> <fct>      <int> <int> <ord>     
## 1 measurement noninoculated inoculated  1.60   StRBOHD        6     6 large     
## 2 measurement noninoculated inoculated -0.0819 StPR1B         6     6 negligible
## 3 measurement noninoculated inoculated -1.44   StCPI8         5     6 large     
## 4 measurement noninoculated inoculated  6.54   StCAB          5     6 large     
## 5 measurement noninoculated inoculated -0.338  StBGLU2        5     6 small     
## 6 measurement noninoculated inoculated -0.0320 StHSP70        6     6 negligible
## 7 measurement noninoculated inoculated -1.97   StPti5         6     5 large     
## 8 measurement noninoculated inoculated -2.36   St13-LOX       6     6 large

##   transcript        group1     group2            p        p.adj p.adj.signif
## 1    StRBOHD noninoculated inoculated 0.0017098192 0.0034196385            *
## 3     StCPI8 noninoculated inoculated 0.0058622374 0.0093795799            *
## 4      StCAB noninoculated inoculated 0.0002442599 0.0006513597            *
## 7     StPti5 noninoculated inoculated 0.0002442599 0.0006513597            *
## 8   St13-LOX noninoculated inoculated 0.0002442599 0.0006513597            *

## ####  ####  
## Distribution tests
## ####  ####  
## # A tibble: 8 × 7
##       n transcript Shapiro_Wilk_BH Anderson_Darling_BH Lilliefors_KS_BH
##   <int> <fct>                <dbl>               <dbl>            <dbl>
## 1    12 StRBOHD           0.0476            0.0339            0.0170   
## 2    12 StPR1B            0.000146          0.00000400        0.0000108
## 3    12 StCPI8            0.0472            0.0339            0.0452   
## 4    10 StCAB             0.0605            0.0777            0.113    
## 5    11 StBGLU2           0.245             0.306             0.440    
## 6    12 StHSP70           0.245             0.377             0.533    
## 7    10 StPti5            0.0221            0.0250            0.0233   
## 8    12 St13-LOX          0.0127            0.00692           0.00302  
## # ℹ 2 more variables: Jarque_Bera_BH <dbl>, DAgostino_Skewness_BH <dbl>
## ####  ####  
## Quantile-Quantile plots
## ####  ####

## ####  ####  
## Test for homogeneity of variance across groups
## ####  ####  
## ####  ####  
## Levene
## ####  ####  
## # A tibble: 8 × 5
##   transcript   df1   df2 statistic       p
##   <fct>      <int> <int>     <dbl>   <dbl>
## 1 StRBOHD        1    10      6.22 0.0318 
## 2 StPR1B         1    10      4.05 0.0720 
## 3 StCPI8         1    10      1.91 0.197  
## 4 StCAB          1     8      9.88 0.0138 
## 5 StBGLU2        1     9     11.0  0.00903
## 6 StHSP70        1    10      2.07 0.181  
## 7 StPti5         1     8      7.63 0.0246 
## 8 St13-LOX       1    10     10.1  0.00975
## ####  ####  
## Brown-Forsythe
## ####  ####  
## # A tibble: 8 × 5
##   transcript   df1   df2 statistic      p
##   <fct>      <int> <int>     <dbl>  <dbl>
## 1 StRBOHD        1    10      4.38 0.0628
## 2 StPR1B         1    10      1.14 0.312 
## 3 StCPI8         1    10      2.30 0.160 
## 4 StCAB          1     8      9.79 0.0140
## 5 StBGLU2        1     9      9.75 0.0123
## 6 StHSP70        1    10      1.42 0.261 
## 7 StPti5         1     8      5.10 0.0539
## 8 St13-LOX       1    10      9.69 0.0110
## ####  ####  
## Fligner
## ####  ####  
## # A tibble: 8 × 5
## # Groups:   transcript [8]
##   transcript statistic p.value parameter method                                 
##   <fct>          <dbl>   <dbl>     <dbl> <chr>                                  
## 1 StRBOHD        3.85  0.0498          1 Fligner-Killeen test of homogeneity of…
## 2 StPR1B         0.633 0.426           1 Fligner-Killeen test of homogeneity of…
## 3 StCPI8         1.79  0.181           1 Fligner-Killeen test of homogeneity of…
## 4 StCAB          5.59  0.0181          1 Fligner-Killeen test of homogeneity of…
## 5 StBGLU2        5.12  0.0236          1 Fligner-Killeen test of homogeneity of…
## 6 StHSP70        0.982 0.322           1 Fligner-Killeen test of homogeneity of…
## 7 StPti5         3.19  0.0739          1 Fligner-Killeen test of homogeneity of…
## 8 St13-LOX       7.59  0.00588         1 Fligner-Killeen test of homogeneity of…
## ####  ####  
## Wilcoxon effect size
## ####  ####  
## # A tibble: 8 × 8
##   .y.         group1        group2     effsize transcript    n1    n2 magnitude
## * <chr>       <chr>         <chr>        <dbl> <fct>      <int> <int> <ord>    
## 1 measurement noninoculated inoculated   0.832 StRBOHD        6     6 large    
## 2 measurement noninoculated inoculated   0.256 StPR1B         6     6 small    
## 3 measurement noninoculated inoculated   0.555 StCPI8         6     6 large    
## 4 measurement noninoculated inoculated   0.809 StCAB          6     4 large    
## 5 measurement noninoculated inoculated   0.220 StBGLU2        6     5 small    
## 6 measurement noninoculated inoculated   0.601 StHSP70        6     6 large    
## 7 measurement noninoculated inoculated   0.826 StPti5         5     5 large    
## 8 measurement noninoculated inoculated   0.832 St13-LOX       6     6 large    
## ####  ####  
## Cohen's d Measure of Effect Size
## ####  ####  
## # A tibble: 8 × 8
##   .y.         group1        group2     effsize transcript    n1    n2 magnitude
## * <chr>       <chr>         <chr>        <dbl> <fct>      <int> <int> <ord>    
## 1 measurement noninoculated inoculated   2.65  StRBOHD        6     6 large    
## 2 measurement noninoculated inoculated   0.702 StPR1B         6     6 moderate 
## 3 measurement noninoculated inoculated   1.31  StCPI8         6     6 large    
## 4 measurement noninoculated inoculated   2.20  StCAB          6     4 large    
## 5 measurement noninoculated inoculated   0.652 StBGLU2        6     5 moderate 
## 6 measurement noninoculated inoculated  -1.55  StHSP70        6     6 large    
## 7 measurement noninoculated inoculated  -3.67  StPti5         5     5 large    
## 8 measurement noninoculated inoculated  -1.88  St13-LOX       6     6 large

##   transcript        group1     group2            p        p.adj p.adj.signif
## 1    StRBOHD noninoculated inoculated 0.0002442599 0.0009770396            *
## 4      StCAB noninoculated inoculated 0.0041524182 0.0083048363            *
## 6    StHSP70 noninoculated inoculated 0.0288226673 0.0461162677            *
## 7     StPti5 noninoculated inoculated 0.0002442599 0.0009770396            *
## 8   St13-LOX noninoculated inoculated 0.0014655594 0.0039081583            *

keep = c("magrittr", "ggplot2", "%!in%", "pal", "my_dir", "fr", "group1_names", "myTable")
all_objs = ls()
keep = c(lsf.str(), intersect(all_objs, keep))
rm(list = setdiff(all_objs, keep))

5 Rywal

5.1 input

Rywal = myTable[myTable$Genotype == 'Rywal', ]

data.table::setDT(Rywal)
Rywal[, SampleID := NULL]
Rywal[, Genotype := NULL]

Rywal$Tissue = as.factor(trimws(Rywal$Tissue))
Rywal$Treatment = factor(trimws(Rywal$Treatment), levels = c("noninoculated", "inoculated"))

5.2 params

groupvars = c("Tissue", "Treatment", "transcript")
measurevar = 'measurement'
dodge = position_dodge(width = 0.5)

myorder = c("StRBOHD", "StPR1B", "StCPI8", "StCAB", "StBGLU2", "StHSP70", "StPti5", "St13-LOX")
group1_names = c("StPti5", "St13-LOX")

5.3 results

cat("Rywal", file = fr, append = TRUE, sep = "\n")
cat("", file = fr, append = TRUE, sep = "\n")


run_analysis_for_strain <- function(data
                                  , strain
                                  , myorder
                                  , pal
                                  , groupvars
                                  , measurevar
                                  , my_dir
                                  , plot_gene_expression
                                  , test_and_plot) {


  strains = trimws(unique(na.omit(data$Strain)))
  temp = data[grep(setdiff(strains, strain), data$Strain, invert = TRUE), ]
  temp[, Strain := NULL]


  shoots = process_data(temp, groupvars, measurevar, scale_treatment = "noninoculated")$shoots
  roots = process_data(temp, groupvars, measurevar, scale_treatment = "noninoculated")$roots
  
  max_per_transcript = shoots %>%
    dplyr::group_by(transcript) %>%
    dplyr::summarise(max_scaled = max(scaled, na.rm = TRUE), .groups = "drop")
  group2_names = setdiff(max_per_transcript$transcript, group1_names)
  y_scales1 = build_y_scales_for(group1_names, max_per_transcript)
  y_scales2 = build_y_scales_for(group2_names, max_per_transcript)
  max_per_transcript = roots %>%
    dplyr::group_by(transcript) %>%
    dplyr::summarise(max_scaled = max(scaled, na.rm = TRUE), .groups = "drop")
  y_scales3 = build_y_scales_for(group1_names, max_per_transcript)
  y_scales4 = build_y_scales_for(group2_names, max_per_transcript)


  

  result_shoots = test_and_plot(data_long_raw = shoots,
                                myorder = myorder,
                                pal = pal,
                                what = paste("shoots", strain),
                                plot_gene_expression_func = plot_gene_expression,
                                groupvars = groupvars,
                                y_scales1 = y_scales1,
                                y_scales2 = y_scales2)
  res = result_shoots$stat.test[, grep("transcript|group1|group2|^p$|p\\.", colnames(result_shoots$stat.test))]
  print(res[res$p.adj.signif != 'ns', ])
  
  cat("", file = fr, append = TRUE, sep = "\n")
  output_text = "permutational t-test"
  cat(output_text, file = fr, append = TRUE, sep = "\n")
  header = paste(colnames(res), collapse = "\t")
  cat(header, file = fr, append = TRUE, sep = "\n")
  output_text = apply(as.data.frame(tibble::as_tibble(res)), 1, function(row) paste(row, collapse = "\t"))
  cat(output_text, file = fr, append = TRUE, sep = "\n")
  cat("", file = fr, append = TRUE, sep = "\n")
  cat("", file = fr, append = TRUE, sep = "\n")
  
  
  ggsave(filename = file.path(my_dir, paste0("Rywal_shoots.", gsub(" ", "", strain), "_1.pdf")),
         plot = result_shoots$plot1, device = pdf, width = 3, height = 8, units = "in", dpi = 900)
  ggsave(filename = file.path(my_dir, paste0("Rywal_shoots.", gsub(" ", "", strain), "_1.svg")),
         plot = result_shoots$plot1, device = svg, width = 3, height = 8, units = "in", dpi = 900)
    
  ggsave(filename = file.path(my_dir, paste0("Rywal_shoots.", gsub(" ", "", strain), "_2.pdf")),
         plot = result_shoots$plot2, device = pdf, width = 9, height = 8, units = "in", dpi = 900)
  ggsave(filename = file.path(my_dir, paste0("Rywal_shoots.", gsub(" ", "", strain), "_2.svg")),
         plot = result_shoots$plot2, device = svg, width = 9, height = 8, units = "in", dpi = 900)
  

  result_roots = test_and_plot(data_long_raw = roots,
                               myorder = myorder,
                               pal = pal,
                               what = paste("roots", strain),
                               plot_gene_expression_func = plot_gene_expression,
                               groupvars = groupvars,
                               y_scales1 = y_scales3,
                               y_scales2 = y_scales4)
  res = result_roots$stat.test[, grep("transcript|group1|group2|^p$|p\\.", colnames(result_roots$stat.test))]
  print(res[res$p.adj.signif != 'ns', ])
  
  cat("", file = fr, append = TRUE, sep = "\n")
  output_text = "permutational t-test"
  cat(output_text, file = fr, append = TRUE, sep = "\n")
  header = paste(colnames(res), collapse = "\t")
  cat(header, file = fr, append = TRUE, sep = "\n")
  output_text = apply(as.data.frame(tibble::as_tibble(res)), 1, function(row) paste(row, collapse = "\t"))
  cat(output_text, file = fr, append = TRUE, sep = "\n")
  cat("", file = fr, append = TRUE, sep = "\n")
  cat("", file = fr, append = TRUE, sep = "\n")
  
  
  ggsave(filename = file.path(my_dir, paste0("Rywal_roots.", gsub(" ", ".", strain), "_1.pdf")),
         plot = result_roots$plot1, device = pdf, width = 3, height = 8, units = "in", dpi = 900)
  ggsave(filename = file.path(my_dir, paste0("Rywal_roots.", gsub(" ", ".", strain), "_1.svg")),
         plot = result_roots$plot1, device = svg, width = 3, height = 8, units = "in", dpi = 900)
  
  ggsave(filename = file.path(my_dir, paste0("Rywal_roots.", gsub(" ", ".", strain), "_2.pdf")),
         plot = result_roots$plot2, device = pdf, width = 9, height = 8, units = "in", dpi = 900)
  ggsave(filename = file.path(my_dir, paste0("Rywal_roots.", gsub(" ", ".", strain), "_2.svg")),
         plot = result_roots$plot2, device = svg, width = 9, height = 8, units = "in", dpi = 900)
  

  list(shoots = result_shoots, roots = result_roots)
}


results_PS216 = run_analysis_for_strain(data = Rywal
                                   , strain = "PS-216"
                                   , myorder
                                   , pal
                                   , groupvars
                                   , measurevar
                                   , my_dir
                                   , plot_gene_expression
                                   , test_and_plot)

## ####  ####  
## Distribution tests
## ####  ####  
## # A tibble: 8 × 7
##       n transcript Shapiro_Wilk_BH Anderson_Darling_BH Lilliefors_KS_BH
##   <int> <fct>                <dbl>               <dbl>            <dbl>
## 1    12 StRBOHD            0.299              0.365             0.724  
## 2    12 StPR1B             0.107              0.0861            0.154  
## 3    12 StCPI8             0.00282            0.000428          0.00514
## 4    10 StCAB              0.107              0.0861            0.154  
## 5    12 StBGLU2            0.107              0.0861            0.154  
## 6    12 StHSP70            0.921              0.715             0.768  
## 7    12 StPti5             0.169              0.268             0.609  
## 8    12 St13-LOX           0.107              0.0861            0.154  
## # ℹ 2 more variables: Jarque_Bera_BH <dbl>, DAgostino_Skewness_BH <dbl>
## ####  ####  
## Quantile-Quantile plots
## ####  ####

## ####  ####  
## Test for homogeneity of variance across groups
## ####  ####  
## ####  ####  
## Levene
## ####  ####  
## # A tibble: 8 × 5
##   transcript   df1   df2 statistic        p
##   <fct>      <int> <int>     <dbl>    <dbl>
## 1 StRBOHD        1    10      2.11 0.177   
## 2 StPR1B         1    10      4.52 0.0595  
## 3 StCPI8         1    10     24.6  0.000567
## 4 StCAB          1     8      3.60 0.0944  
## 5 StBGLU2        1    10      5.34 0.0435  
## 6 StHSP70        1    10      3.15 0.106   
## 7 StPti5         1    10      2.91 0.119   
## 8 St13-LOX       1    10      4.12 0.0698  
## ####  ####  
## Brown-Forsythe
## ####  ####  
## # A tibble: 8 × 5
##   transcript   df1   df2 statistic      p
##   <fct>      <int> <int>     <dbl>  <dbl>
## 1 StRBOHD        1    10     1.86  0.202 
## 2 StPR1B         1    10     3.35  0.0970
## 3 StCPI8         1    10     6.58  0.0281
## 4 StCAB          1     8     3.08  0.117 
## 5 StBGLU2        1    10     0.842 0.380 
## 6 StHSP70        1    10     3.06  0.111 
## 7 StPti5         1    10     0.494 0.498 
## 8 St13-LOX       1    10     2.13  0.175 
## ####  ####  
## Fligner
## ####  ####  
## # A tibble: 8 × 5
## # Groups:   transcript [8]
##   transcript statistic p.value parameter method                                 
##   <fct>          <dbl>   <dbl>     <dbl> <chr>                                  
## 1 StRBOHD      1.42    0.234           1 Fligner-Killeen test of homogeneity of…
## 2 StPR1B       3.41    0.0649          1 Fligner-Killeen test of homogeneity of…
## 3 StCPI8       7.57    0.00592         1 Fligner-Killeen test of homogeneity of…
## 4 StCAB        3.01    0.0826          1 Fligner-Killeen test of homogeneity of…
## 5 StBGLU2      0.254   0.614           1 Fligner-Killeen test of homogeneity of…
## 6 StHSP70      3.38    0.0660          1 Fligner-Killeen test of homogeneity of…
## 7 StPti5       0.00225 0.962           1 Fligner-Killeen test of homogeneity of…
## 8 St13-LOX     3.57    0.0588          1 Fligner-Killeen test of homogeneity of…
## ####  ####  
## Wilcoxon effect size
## ####  ####  
## # A tibble: 8 × 8
##   .y.         group1        group2     effsize transcript    n1    n2 magnitude
## * <chr>       <chr>         <chr>        <dbl> <fct>      <int> <int> <ord>    
## 1 measurement noninoculated inoculated  0.832  StRBOHD        6     6 large    
## 2 measurement noninoculated inoculated  0.0462 StPR1B         6     6 small    
## 3 measurement noninoculated inoculated  0.832  StCPI8         6     6 large    
## 4 measurement noninoculated inoculated  0.809  StCAB          4     6 large    
## 5 measurement noninoculated inoculated  0      StBGLU2        6     6 small    
## 6 measurement noninoculated inoculated  0.277  StHSP70        6     6 small    
## 7 measurement noninoculated inoculated  0.324  StPti5         6     6 moderate 
## 8 measurement noninoculated inoculated  0.832  St13-LOX       6     6 large    
## ####  ####  
## Cohen's d Measure of Effect Size
## ####  ####  
## # A tibble: 8 × 8
##   .y.         group1        group2     effsize transcript    n1    n2 magnitude
## * <chr>       <chr>         <chr>        <dbl> <fct>      <int> <int> <ord>    
## 1 measurement noninoculated inoculated   4.19  StRBOHD        6     6 large    
## 2 measurement noninoculated inoculated   0.362 StPR1B         6     6 small    
## 3 measurement noninoculated inoculated  -1.45  StCPI8         6     6 large    
## 4 measurement noninoculated inoculated   3.13  StCAB          4     6 large    
## 5 measurement noninoculated inoculated  -0.215 StBGLU2        6     6 small    
## 6 measurement noninoculated inoculated   0.556 StHSP70        6     6 moderate 
## 7 measurement noninoculated inoculated  -0.209 StPti5         6     6 small    
## 8 measurement noninoculated inoculated  -4.28  St13-LOX       6     6 large

##   transcript        group1     group2            p        p.adj p.adj.signif
## 1    StRBOHD noninoculated inoculated 0.0002442599 0.0004885198            *
## 3     StCPI8 noninoculated inoculated 0.0002442599 0.0004885198            *
## 4      StCAB noninoculated inoculated 0.0002442599 0.0004885198            *
## 8   St13-LOX noninoculated inoculated 0.0002442599 0.0004885198            *

## ####  ####  
## Distribution tests
## ####  ####  
## # A tibble: 8 × 7
##       n transcript Shapiro_Wilk_BH Anderson_Darling_BH Lilliefors_KS_BH
##   <int> <fct>                <dbl>               <dbl>            <dbl>
## 1    11 StRBOHD           0.324             0.332              0.330   
## 2    12 StPR1B            0.000407          0.00000455         0.000173
## 3    10 StCPI8            0.164             0.249              0.337   
## 4     9 StCAB             0.0306            0.0390             0.171   
## 5    12 StBGLU2           0.00109           0.000552           0.00929 
## 6    12 StHSP70           0.155             0.126              0.271   
## 7    12 StPti5            0.00284           0.000925           0.00929 
## 8    12 St13-LOX          0.0312            0.0474             0.0632  
## # ℹ 2 more variables: Jarque_Bera_BH <dbl>, DAgostino_Skewness_BH <dbl>
## ####  ####  
## Quantile-Quantile plots
## ####  ####

## ####  ####  
## Test for homogeneity of variance across groups
## ####  ####  
## ####  ####  
## Levene
## ####  ####  
## # A tibble: 8 × 5
##   transcript   df1   df2 statistic        p
##   <fct>      <int> <int>     <dbl>    <dbl>
## 1 StRBOHD        1     9     4.57  0.0613  
## 2 StPR1B         1    10    29.5   0.000290
## 3 StCPI8         1     8     0.779 0.403   
## 4 StCAB          1     7     0.148 0.712   
## 5 StBGLU2        1    10     4.74  0.0546  
## 6 StHSP70        1    10     0.396 0.543   
## 7 StPti5         1    10    15.7   0.00268 
## 8 St13-LOX       1    10     5.68  0.0384  
## ####  ####  
## Brown-Forsythe
## ####  ####  
## # A tibble: 8 × 5
##   transcript   df1   df2 statistic       p
##   <fct>      <int> <int>     <dbl>   <dbl>
## 1 StRBOHD        1     9    4.40   0.0653 
## 2 StPR1B         1    10    4.59   0.0579 
## 3 StCPI8         1     8    0.0906 0.771  
## 4 StCAB          1     7    0.231  0.645  
## 5 StBGLU2        1    10    3.81   0.0795 
## 6 StHSP70        1    10    0.0184 0.895  
## 7 StPti5         1    10   10.7    0.00851
## 8 St13-LOX       1    10    4.40   0.0623 
## ####  ####  
## Fligner
## ####  ####  
## # A tibble: 8 × 5
## # Groups:   transcript [8]
##   transcript statistic p.value parameter method                                 
##   <fct>          <dbl>   <dbl>     <dbl> <chr>                                  
## 1 StRBOHD       3.21   0.0734          1 Fligner-Killeen test of homogeneity of…
## 2 StPR1B        7.61   0.00579         1 Fligner-Killeen test of homogeneity of…
## 3 StCPI8        0.0695 0.792           1 Fligner-Killeen test of homogeneity of…
## 4 StCAB         0.466  0.495           1 Fligner-Killeen test of homogeneity of…
## 5 StBGLU2       7.59   0.00588         1 Fligner-Killeen test of homogeneity of…
## 6 StHSP70       0.0916 0.762           1 Fligner-Killeen test of homogeneity of…
## 7 StPti5        7.59   0.00588         1 Fligner-Killeen test of homogeneity of…
## 8 St13-LOX      7.57   0.00592         1 Fligner-Killeen test of homogeneity of…
## ####  ####  
## Wilcoxon effect size
## ####  ####  
## # A tibble: 8 × 8
##   .y.         group1        group2     effsize transcript    n1    n2 magnitude
## * <chr>       <chr>         <chr>        <dbl> <fct>      <int> <int> <ord>    
## 1 measurement noninoculated inoculated   0.661 StRBOHD        6     5 large    
## 2 measurement noninoculated inoculated   0.603 StPR1B         6     6 large    
## 3 measurement noninoculated inoculated   0.495 StCPI8         5     5 moderate 
## 4 measurement noninoculated inoculated   0.245 StCAB          5     4 small    
## 5 measurement noninoculated inoculated   0.601 StBGLU2        6     6 large    
## 6 measurement noninoculated inoculated   0.324 StHSP70        6     6 moderate 
## 7 measurement noninoculated inoculated   0.832 StPti5         6     6 large    
## 8 measurement noninoculated inoculated   0.832 St13-LOX       6     6 large    
## ####  ####  
## Cohen's d Measure of Effect Size
## ####  ####  
## # A tibble: 8 × 8
##   .y.         group1        group2     effsize transcript    n1    n2 magnitude
## * <chr>       <chr>         <chr>        <dbl> <fct>      <int> <int> <ord>    
## 1 measurement noninoculated inoculated   1.70  StRBOHD        6     5 large    
## 2 measurement noninoculated inoculated  -1.08  StPR1B         6     6 large    
## 3 measurement noninoculated inoculated  -1.37  StCPI8         5     5 large    
## 4 measurement noninoculated inoculated  -0.254 StCAB          5     4 small    
## 5 measurement noninoculated inoculated  -1.09  StBGLU2        6     6 large    
## 6 measurement noninoculated inoculated  -0.712 StHSP70        6     6 moderate 
## 7 measurement noninoculated inoculated  -1.60  StPti5         6     6 large    
## 8 measurement noninoculated inoculated  -2.41  St13-LOX       6     6 large

##   transcript        group1     group2            p       p.adj p.adj.signif
## 1    StRBOHD noninoculated inoculated 0.0236932096 0.037909135            *
## 2     StPR1B noninoculated inoculated 0.0183194919 0.037909135            *
## 5    StBGLU2 noninoculated inoculated 0.0195407914 0.037909135            *
## 7     StPti5 noninoculated inoculated 0.0002442599 0.001954079            *
## 8   St13-LOX noninoculated inoculated 0.0019540791 0.007816317            *

results_PS68 = run_analysis_for_strain(data = Rywal
                                   , strain = "PS-68"
                                   , myorder
                                   , pal
                                   , groupvars
                                   , measurevar
                                   , my_dir
                                   , plot_gene_expression
                                   , test_and_plot)

## ####  ####  
## Distribution tests
## ####  ####  
## # A tibble: 8 × 7
##       n transcript Shapiro_Wilk_BH Anderson_Darling_BH Lilliefors_KS_BH
##   <int> <fct>                <dbl>               <dbl>            <dbl>
## 1    12 StRBOHD            0.538               0.519            0.424  
## 2    12 StPR1B             0.00485             0.00198          0.0133 
## 3    12 StCPI8             0.00485             0.00198          0.00771
## 4    10 StCAB              0.0381              0.0405           0.173  
## 5    12 StBGLU2            0.0135              0.0307           0.129  
## 6    12 StHSP70            0.538               0.519            0.250  
## 7    11 StPti5             0.0225              0.0307           0.132  
## 8    12 St13-LOX           0.116               0.111            0.173  
## # ℹ 2 more variables: Jarque_Bera_BH <dbl>, DAgostino_Skewness_BH <dbl>
## ####  ####  
## Quantile-Quantile plots
## ####  ####

## ####  ####  
## Test for homogeneity of variance across groups
## ####  ####  
## ####  ####  
## Levene
## ####  ####  
## # A tibble: 8 × 5
##   transcript   df1   df2 statistic        p
##   <fct>      <int> <int>     <dbl>    <dbl>
## 1 StRBOHD        1    10    0.0204 0.889   
## 2 StPR1B         1    10   14.2    0.00364 
## 3 StCPI8         1    10   27.2    0.000393
## 4 StCAB          1     8    4.46   0.0676  
## 5 StBGLU2        1    10    2.46   0.148   
## 6 StHSP70        1    10   10.6    0.00863 
## 7 StPti5         1     9    0.104  0.755   
## 8 St13-LOX       1    10    3.34   0.0974  
## ####  ####  
## Brown-Forsythe
## ####  ####  
## # A tibble: 8 × 5
##   transcript   df1   df2 statistic        p
##   <fct>      <int> <int>     <dbl>    <dbl>
## 1 StRBOHD        1    10    0.0150 0.905   
## 2 StPR1B         1    10    2.64   0.135   
## 3 StCPI8         1    10   24.8    0.000552
## 4 StCAB          1     8    4.02   0.0800  
## 5 StBGLU2        1    10    1.62   0.232   
## 6 StHSP70        1    10   10.1    0.00984 
## 7 StPti5         1     9    0.0408 0.844   
## 8 St13-LOX       1    10    2.74   0.129   
## ####  ####  
## Fligner
## ####  ####  
## # A tibble: 8 × 5
## # Groups:   transcript [8]
##   transcript statistic p.value parameter method                                 
##   <fct>          <dbl>   <dbl>     <dbl> <chr>                                  
## 1 StRBOHD       0.0720 0.788           1 Fligner-Killeen test of homogeneity of…
## 2 StPR1B        2.86   0.0909          1 Fligner-Killeen test of homogeneity of…
## 3 StCPI8        7.57   0.00592         1 Fligner-Killeen test of homogeneity of…
## 4 StCAB         4.52   0.0336          1 Fligner-Killeen test of homogeneity of…
## 5 StBGLU2       1.79   0.180           1 Fligner-Killeen test of homogeneity of…
## 6 StHSP70       6.09   0.0136          1 Fligner-Killeen test of homogeneity of…
## 7 StPti5        0.355  0.551           1 Fligner-Killeen test of homogeneity of…
## 8 St13-LOX      4.16   0.0415          1 Fligner-Killeen test of homogeneity of…
## ####  ####  
## Wilcoxon effect size
## ####  ####  
## # A tibble: 8 × 8
##   .y.         group1        group2     effsize transcript    n1    n2 magnitude
## * <chr>       <chr>         <chr>        <dbl> <fct>      <int> <int> <ord>    
## 1 measurement noninoculated inoculated   0.786 StRBOHD        6     6 large    
## 2 measurement noninoculated inoculated   0.277 StPR1B         6     6 small    
## 3 measurement noninoculated inoculated   0.693 StCPI8         6     6 large    
## 4 measurement noninoculated inoculated   0.809 StCAB          4     6 large    
## 5 measurement noninoculated inoculated   0.370 StBGLU2        6     6 moderate 
## 6 measurement noninoculated inoculated   0.231 StHSP70        6     6 small    
## 7 measurement noninoculated inoculated   0.385 StPti5         6     5 moderate 
## 8 measurement noninoculated inoculated   0.832 St13-LOX       6     6 large    
## ####  ####  
## Cohen's d Measure of Effect Size
## ####  ####  
## # A tibble: 8 × 8
##   .y.         group1        group2     effsize transcript    n1    n2 magnitude
## * <chr>       <chr>         <chr>        <dbl> <fct>      <int> <int> <ord>    
## 1 measurement noninoculated inoculated   2.81  StRBOHD        6     6 large    
## 2 measurement noninoculated inoculated  -0.862 StPR1B         6     6 large    
## 3 measurement noninoculated inoculated  -1.47  StCPI8         6     6 large    
## 4 measurement noninoculated inoculated   3.32  StCAB          4     6 large    
## 5 measurement noninoculated inoculated  -0.722 StBGLU2        6     6 moderate 
## 6 measurement noninoculated inoculated  -0.489 StHSP70        6     6 small    
## 7 measurement noninoculated inoculated  -0.497 StPti5         6     5 small    
## 8 measurement noninoculated inoculated  -3.64  St13-LOX       6     6 large

##   transcript        group1     group2            p       p.adj p.adj.signif
## 1    StRBOHD noninoculated inoculated 0.0014655594 0.005862237            *
## 3     StCPI8 noninoculated inoculated 0.0092818759 0.018563752            *
## 4      StCAB noninoculated inoculated 0.0043966781 0.011724475            *
## 8   St13-LOX noninoculated inoculated 0.0002442599 0.001954079            *

## ####  ####  
## Distribution tests
## ####  ####  
## # A tibble: 8 × 7
##       n transcript Shapiro_Wilk_BH Anderson_Darling_BH Lilliefors_KS_BH
##   <int> <fct>                <dbl>               <dbl>            <dbl>
## 1    12 StRBOHD           0.605              0.741           0.878     
## 2    12 StPR1B            0.000364           0.0000123       0.00000862
## 3    11 StCPI8            0.0915             0.166           0.432     
## 4     9 StCAB             0.199              0.240           0.244     
## 5    12 StBGLU2           0.00166            0.000675        0.000800  
## 6    12 StHSP70           0.0435             0.0793          0.244     
## 7    12 StPti5            0.0435             0.0424          0.0148    
## 8    12 St13-LOX          0.0184             0.0105          0.0148    
## # ℹ 2 more variables: Jarque_Bera_BH <dbl>, DAgostino_Skewness_BH <dbl>
## ####  ####  
## Quantile-Quantile plots
## ####  ####

## ####  ####  
## Test for homogeneity of variance across groups
## ####  ####  
## ####  ####  
## Levene
## ####  ####  
## # A tibble: 8 × 5
##   transcript   df1   df2 statistic       p
##   <fct>      <int> <int>     <dbl>   <dbl>
## 1 StRBOHD        1    10     1.95  0.193  
## 2 StPR1B         1    10     7.75  0.0193 
## 3 StCPI8         1     9     2.79  0.129  
## 4 StCAB          1     7     0.758 0.413  
## 5 StBGLU2        1    10     7.26  0.0225 
## 6 StHSP70        1    10     0.343 0.571  
## 7 StPti5         1    10    10.6   0.00866
## 8 St13-LOX       1    10    16.3   0.00236
## ####  ####  
## Brown-Forsythe
## ####  ####  
## # A tibble: 8 × 5
##   transcript   df1   df2 statistic       p
##   <fct>      <int> <int>     <dbl>   <dbl>
## 1 StRBOHD        1    10     1.93  0.195  
## 2 StPR1B         1    10     5.86  0.0360 
## 3 StCPI8         1     9     2.64  0.139  
## 4 StCAB          1     7     0.247 0.634  
## 5 StBGLU2        1    10     4.06  0.0715 
## 6 StHSP70        1    10     0.118 0.738  
## 7 StPti5         1    10     6.73  0.0268 
## 8 St13-LOX       1    10    14.5   0.00347
## ####  ####  
## Fligner
## ####  ####  
## # A tibble: 8 × 5
## # Groups:   transcript [8]
##   transcript statistic p.value parameter method                                 
##   <fct>          <dbl>   <dbl>     <dbl> <chr>                                  
## 1 StRBOHD       1.63   0.202           1 Fligner-Killeen test of homogeneity of…
## 2 StPR1B        7.61   0.00579         1 Fligner-Killeen test of homogeneity of…
## 3 StCPI8        1.38   0.241           1 Fligner-Killeen test of homogeneity of…
## 4 StCAB         0.0363 0.849           1 Fligner-Killeen test of homogeneity of…
## 5 StBGLU2       7.57   0.00592         1 Fligner-Killeen test of homogeneity of…
## 6 StHSP70       0.0202 0.887           1 Fligner-Killeen test of homogeneity of…
## 7 StPti5        7.59   0.00588         1 Fligner-Killeen test of homogeneity of…
## 8 St13-LOX      6.09   0.0136          1 Fligner-Killeen test of homogeneity of…
## ####  ####  
## Wilcoxon effect size
## ####  ####  
## # A tibble: 8 × 8
##   .y.         group1        group2     effsize transcript    n1    n2 magnitude
## * <chr>       <chr>         <chr>        <dbl> <fct>      <int> <int> <ord>    
## 1 measurement noninoculated inoculated   0.508 StRBOHD        6     6 large    
## 2 measurement noninoculated inoculated   0.488 StPR1B         6     6 moderate 
## 3 measurement noninoculated inoculated   0.661 StCPI8         5     6 large    
## 4 measurement noninoculated inoculated   0.327 StCAB          5     4 moderate 
## 5 measurement noninoculated inoculated   0.370 StBGLU2        6     6 moderate 
## 6 measurement noninoculated inoculated   0.139 StHSP70        6     6 small    
## 7 measurement noninoculated inoculated   0.832 StPti5         6     6 large    
## 8 measurement noninoculated inoculated   0.832 St13-LOX       6     6 large    
## ####  ####  
## Cohen's d Measure of Effect Size
## ####  ####  
## # A tibble: 8 × 8
##   .y.         group1        group2     effsize transcript    n1    n2 magnitude
## * <chr>       <chr>         <chr>        <dbl> <fct>      <int> <int> <ord>    
## 1 measurement noninoculated inoculated   1.08  StRBOHD        6     6 large    
## 2 measurement noninoculated inoculated  -1.01  StPR1B         6     6 large    
## 3 measurement noninoculated inoculated  -1.68  StCPI8         5     6 large    
## 4 measurement noninoculated inoculated  -0.225 StCAB          5     4 small    
## 5 measurement noninoculated inoculated  -0.952 StBGLU2        6     6 large    
## 6 measurement noninoculated inoculated  -0.340 StHSP70        6     6 small    
## 7 measurement noninoculated inoculated  -3.44  StPti5         6     6 large    
## 8 measurement noninoculated inoculated  -9.33  St13-LOX       6     6 large

##   transcript        group1     group2            p        p.adj p.adj.signif
## 7     StPti5 noninoculated inoculated 0.0002442599 0.0009770396            *
## 8   St13-LOX noninoculated inoculated 0.0002442599 0.0009770396            *

6 sessionInfo

sessionInfo()
## R version 4.4.1 (2024-06-14 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 11 x64 (build 26100)
## 
## Matrix products: default
## 
## 
## locale:
## [1] LC_COLLATE=English_United Kingdom.utf8 
## [2] LC_CTYPE=English_United Kingdom.utf8   
## [3] LC_MONETARY=English_United Kingdom.utf8
## [4] LC_NUMERIC=C                           
## [5] LC_TIME=English_United Kingdom.utf8    
## 
## time zone: Europe/Ljubljana
## tzcode source: internal
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] ggplot2_4.0.0  magrittr_2.0.3
## 
## loaded via a namespace (and not attached):
##   [1] DBI_1.2.3             Rdpack_2.6.4          sandwich_3.1-1       
##   [4] MKinfer_1.2           rlang_1.1.6           multcomp_1.4-28      
##   [7] miceadds_3.18-36      tseries_0.10-58       matrixStats_1.5.0    
##  [10] compiler_4.4.1        vctrs_0.6.5           quadprog_1.5-8       
##  [13] pkgconfig_2.0.3       shape_1.4.6.1         crayon_1.5.3         
##  [16] fastmap_1.2.0         backports_1.5.0       labeling_0.4.3       
##  [19] utf8_1.2.6            rmarkdown_2.30        exactRankTests_0.8-35
##  [22] nloptr_2.2.1          purrr_1.0.4           xfun_0.53            
##  [25] glmnet_4.1-10         modeltools_0.2-24     jomo_2.7-6           
##  [28] cachem_1.1.0          jsonlite_2.0.0        gmp_0.7-5            
##  [31] pan_1.9               broom_1.0.10          parallel_4.4.1       
##  [34] R6_2.6.1              coin_1.4-3            bslib_0.9.0          
##  [37] stringi_1.8.7         RColorBrewer_1.1-3    rpart_4.1.24         
##  [40] car_3.1-3             boot_1.3-31           jquerylib_0.1.4      
##  [43] Rcpp_1.0.14           iterators_1.0.14      knitr_1.50           
##  [46] zoo_1.8-14            nnet_7.3-20           Matrix_1.7-1         
##  [49] splines_4.4.1         tidyselect_1.2.1      rstudioapi_0.17.1    
##  [52] dichromat_2.0-0.1     abind_1.4-8           yaml_2.3.10          
##  [55] codetools_0.2-20      curl_6.2.2            arrangements_1.1.9   
##  [58] lattice_0.22-6        tibble_3.3.0          quantmod_0.4.28      
##  [61] withr_3.0.2           S7_0.2.0              evaluate_1.0.5       
##  [64] moments_0.14.1        survival_3.8-3        zip_2.3.3            
##  [67] xts_0.14.1            pillar_1.11.1         ggpubr_0.6.1         
##  [70] carData_3.0-5         mice_3.18.0           nortest_1.0-4        
##  [73] foreach_1.5.2         stats4_4.4.1          reformulas_0.4.1     
##  [76] generics_0.1.4        TTR_0.24.4            scales_1.4.0         
##  [79] minqa_1.2.8           glue_1.8.0            tools_4.4.1          
##  [82] data.table_1.17.8     MKdescr_0.9           lme4_1.1-37          
##  [85] openxlsx_4.2.8        ggsignif_0.6.4        mvtnorm_1.3-3        
##  [88] cowplot_1.2.0         grid_4.4.1            mitools_2.4          
##  [91] tidyr_1.3.1           rbibutils_2.3         libcoin_1.0-10       
##  [94] nlme_3.1-166          Formula_1.2-5         cli_3.6.5            
##  [97] dplyr_1.1.4           ggh4x_0.3.1           gtable_0.3.6         
## [100] rstatix_0.7.2         sass_0.4.10           digest_0.6.37        
## [103] TH.data_1.1-4         farver_2.1.2          htmltools_0.5.8.1    
## [106] lifecycle_1.0.4       mitml_0.4-5           MASS_7.3-64